. Set working directory. Load packages. Set theme for figures

setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) ## Data must be in same folder

if (!require(librarian)) {
  install.packages("librarian")
  library(librarian)
}
## Loading required package: librarian
shelf(tidyverse, ggplot2, ggpmisc, viridis, ggpubr, lme4, lmerTest, DHARMa, sjPlot, glmmTMB, effects, ggeffects, AICcmodavg, modelsummary, emmeans, data.table, MuMIn, Durga, car, truncnorm, random, performance, png, grid, loo, bayesplot, brms, png, grid, gtools, ggnewscale)
## 
##   The 'cran_repo' argument in shelf() was not set, so it will use
##   cran_repo = 'https://cran.r-project.org' by default.
## 
##   To avoid this message, set the 'cran_repo' argument to a CRAN
##   mirror URL (see https://cran.r-project.org/mirrors.html) or set
##   'quiet = TRUE'.
## Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch: 
## glmmTMB was built with TMB package version 1.9.15
## Current TMB package version is 1.9.16
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
## Warning: package 'bayesplot' was built under R version 4.4.3
## Warning: package 'ggnewscale' was built under R version 4.4.3
set_theme(base = theme_bw(),
          axis.title.size = 1.5,
          axis.textsize = 1.0, axis.title.color = 'black', axis.textcolor = 'black',
          legend.title.size = 1.5)

. Load datasets

df = dataset on chain size and force allchainposdf = ant posture dataset thoraxDF = thorax length dataset

df <- read.csv("ChainData.csv") 
allchainposdf <- read.csv("ChainPos.csv") 
thoraxDF <- read.csv("ThoraxLengths.csv")

. Summarise dataset by arrangement and calculate chain length ratio

Arrangement changes every time an ant joins or leaves the pulling effort. Note that ‘Chain length ratio’ is named ‘propmax’ throughout the code.

meandf <- df %>%
  group_by(replicate, newarr) %>%
  mutate(n_row = n()) %>% ## Duration of arrangement in seconds
  mutate(
    state = ifelse(ptime < 0, 'grow', 'decay'),
    meanforceN = mean(forceN),
    meanforceperantN = mean(forceperantN)) %>%
  slice_head() %>%
  mutate(
    prop1 = chain1*1/total.ants,
    prop2 = chain2*2/total.ants,
    prop3 = chain3*3/total.ants,
    prop4 = chain4*4/total.ants,
    prop5 = chain5*5/total.ants,
  ) %>%
  mutate(
    propmax = as.factor(which.max(c(prop1, prop2, prop3, prop4, prop5))), ## chain length ratio
    state = factor(state, levels = c("grow","decay"))
  )

. Summary of experimental observations

length(unique(meandf$replicate)) ## 15 reps
## [1] 15
length(unique(meandf$arrangement)) ## 108 arrangements
## [1] 108
length(unique(meandf$colony)) ## 5 colonies
## [1] 5
n <- meandf %>% group_by(replicate) %>% summarise(n = length(unique(newarr))); sum(n$n) ## 364 new arrivals
## [1] 364
nrow(meandf) ## 364 observations in total
## [1] 364

. Average force produced by ants

avg_forceperant <- mean(meandf$meanforceperantN) ## Average individual force output throughout all replicates

sem_forceperant <- sd(meandf$meanforceperantN)/sqrt(length(meandf$meanforceperantN)) ## Standard deviation

max_forceperant <- max(meandf$meanforceperantN) ## Maximum individual force output recorded

mean(meandf[meandf$total.ants == 1,]$meanforceperantN) ## Mean individual force output when team size == 1
## [1] 2.891822
mean(meandf[meandf$total.ants == 15,]$meanforceperantN) ## Mean individual force output when team size == 15
## [1] 5.064689

. Data analysis and visualization

. GLMM on team force output

# Fit model
fit1 <- glmmTMB(sqrt(meanforceN) ~ total.ants * state + propmax + (1|colony/replicate), data = meandf) 

# Diagnostics
check_model(fit1)
## `check_outliers()` does not yet support models of class `glmmTMB`.

# Results
summary(fit1)
##  Family: gaussian  ( identity )
## Formula:          
## sqrt(meanforceN) ~ total.ants * state + propmax + (1 | colony/replicate)
## Data: meandf
## 
##      AIC      BIC   logLik deviance df.resid 
##   1065.5   1104.5   -522.8   1045.5      354 
## 
## Random effects:
## 
## Conditional model:
##  Groups           Name        Variance Std.Dev.
##  replicate:colony (Intercept) 0.4584   0.6770  
##  colony           (Intercept) 0.2183   0.4673  
##  Residual                     0.9259   0.9623  
## Number of obs: 364, groups:  replicate:colony, 15; colony, 5
## 
## Dispersion estimate for gaussian family (sigma^2): 0.926 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.006404   0.350017   0.018 0.985403    
## total.ants             0.586778   0.023665  24.795  < 2e-16 ***
## statedecay             1.107256   0.251384   4.405 1.06e-05 ***
## propmax2               0.463790   0.160862   2.883 0.003937 ** 
## propmax3               0.601826   0.180135   3.341 0.000835 ***
## propmax4               0.583375   0.322766   1.807 0.070696 .  
## total.ants:statedecay -0.045518   0.027950  -1.629 0.103406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(fit1, "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: sqrt(meanforceN)
##                     Chisq Df Pr(>Chisq)    
## (Intercept)        0.0003  1    0.98540    
## total.ants       614.7799  1  < 2.2e-16 ***
## state             19.4009  1   1.06e-05 ***
## propmax           13.0406  3    0.00455 ** 
## total.ants:state   2.6522  1    0.10341    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Forest plot
plot_model(fit1)

## Pairwise comparison between chain length ratio 
emms <- emmeans(fit1, ~ propmax)
pairs(emms)
##  contrast            estimate    SE  df t.ratio p.value
##  propmax1 - propmax2  -0.4638 0.161 354  -2.883  0.0216
##  propmax1 - propmax3  -0.6018 0.180 354  -3.341  0.0051
##  propmax1 - propmax4  -0.5834 0.323 354  -1.807  0.2715
##  propmax2 - propmax3  -0.1380 0.154 354  -0.894  0.8081
##  propmax2 - propmax4  -0.1196 0.324 354  -0.369  0.9828
##  propmax3 - propmax4   0.0185 0.322 354   0.057  0.9999
## 
## Results are averaged over the levels of: state 
## Note: contrasts are still on the sqrt scale. Consider using
##       regrid() if you want contrasts of back-transformed estimates. 
## P value adjustment: tukey method for comparing a family of 4 estimates
### PLOTTING ###

## Predict based on total number of ants and state
dat <- ggpredict(fit1, terms = c("total.ants [all]", "state [all]")) %>%
  rename(state = group)
## Model has sqrt-transformed response. Back-transforming predictions to
##   original response scale. Standard errors are still on the transformed
##   scale.
## Visualise results 
FigX1 <- ggplot() +
  geom_point(data = meandf, aes(total.ants, meanforceN, color = state, fill = state), size = 3) +
  geom_line(data = dat, aes(x, predicted, color = state)) +
  geom_ribbon(data = dat, aes(x = x, ymin = conf.low, ymax = conf.high, fill = state), alpha = .4) +
  labs(
    x = "Total nº of ants",
    y = "Total Force (mN)",
    fill = "Phase",
    colour = "Phase") +
  scale_color_manual(values = c("#00BFC4", "#F8766D"), labels = c("Growth", "Decay")) +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"), labels = c("Growth", "Decay")) +
  theme_pubr(legend = c(0.1, 0.9),
             base_size = 13); FigX1

## Effect without transformation
dat1 <- ggeffect(fit1, terms = c("total.ants [all]")) %>%
  rename(state = group) 

ggplot(meandf) +
  geom_point(aes(total.ants, sqrt(meanforceN), color = state)) +
  geom_line(data = dat1, aes(x, predicted)) +
  geom_ribbon(data = dat1, aes(x = x, ymin = conf.low, ymax = conf.high), alpha = .2) +
  labs(
    x = "Total nº of ants",
    y = "Square-root Mean Force (mN)"
  )

## Durga plot
d <- Durga::DurgaDiff(meandf, "meanforceN", "propmax")
DurgaPlot(d, contrasts = c("2 - 1", "3 - 2", "4 - 3"), xlab = "PropMax",left.ylab = "Mean force (mN)") 

. GLMM on individual force output

# Fit model
fit2 <- glmmTMB(sqrt(meanforceperantN) ~ total.ants * state + propmax + (1|colony/replicate), data = meandf) 

# Diagnostics
check_model(fit2)
## `check_outliers()` does not yet support models of class `glmmTMB`.

# Results
summary(fit2)
##  Family: gaussian  ( identity )
## Formula:          sqrt(meanforceperantN) ~ total.ants * state + propmax + (1 |  
##     colony/replicate)
## Data: meandf
## 
##      AIC      BIC   logLik deviance df.resid 
##    361.1    400.0   -170.5    341.1      354 
## 
## Random effects:
## 
## Conditional model:
##  Groups           Name        Variance Std.Dev.
##  replicate:colony (Intercept) 0.0702   0.2650  
##  colony           (Intercept) 0.0313   0.1769  
##  Residual                     0.1334   0.3653  
## Number of obs: 364, groups:  replicate:colony, 15; colony, 5
## 
## Dispersion estimate for gaussian family (sigma^2): 0.133 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.871253   0.133909   6.506 7.70e-11 ***
## total.ants             0.096704   0.009009  10.734  < 2e-16 ***
## statedecay             0.649111   0.095508   6.796 1.07e-11 ***
## propmax2               0.144668   0.061104   2.368 0.017906 *  
## propmax3               0.220739   0.068415   3.226 0.001253 ** 
## propmax4               0.199508   0.122568   1.628 0.103582    
## total.ants:statedecay -0.040088   0.010615  -3.777 0.000159 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(fit2, "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: sqrt(meanforceperantN)
##                    Chisq Df Pr(>Chisq)    
## (Intercept)       42.332  1  7.701e-11 ***
## total.ants       115.222  1  < 2.2e-16 ***
## state             46.191  1  1.072e-11 ***
## propmax           11.202  3   0.010681 *  
## total.ants:state  14.263  1   0.000159 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Forest plot
plot_model(fit2)

# Pairwise comparison between chain length ratio 
emms.ifo <- emmeans(fit2, ~ propmax)
pairs(emms.ifo)
##  contrast            estimate     SE  df t.ratio p.value
##  propmax1 - propmax2  -0.1447 0.0611 354  -2.368  0.0853
##  propmax1 - propmax3  -0.2207 0.0684 354  -3.226  0.0075
##  propmax1 - propmax4  -0.1995 0.1230 354  -1.628  0.3643
##  propmax2 - propmax3  -0.0761 0.0586 354  -1.299  0.5642
##  propmax2 - propmax4  -0.0548 0.1230 354  -0.445  0.9705
##  propmax3 - propmax4   0.0212 0.1220 354   0.174  0.9981
## 
## Results are averaged over the levels of: state 
## Note: contrasts are still on the sqrt scale. Consider using
##       regrid() if you want contrasts of back-transformed estimates. 
## P value adjustment: tukey method for comparing a family of 4 estimates
### PLOTTING ###
dat2 <- ggpredict(fit2, terms = c("total.ants [all]", "state [all]")) %>%
  rename(state = group)
## Model has sqrt-transformed response. Back-transforming predictions to
##   original response scale. Standard errors are still on the transformed
##   scale.
## Effect of total number of ants and state
FigX2 <- ggplot() +
  geom_jitter(data = meandf, aes(total.ants, meanforceperantN, color = state), height = 0, width = 0) +
  scale_x_continuous(breaks = c(1, 5, 10, 15)) +
  geom_line(data = dat2, aes(x, predicted, color = state)) +
  geom_ribbon(data = dat2, aes(x = x, ymin = conf.low, ymax = conf.high, fill = state), alpha = .2) +
  labs(
    x = "Total nº of ants",
    y = "Mean Force per ant (mN)",
    colour = "Phase",
    fill = "Phase"
  ) +
  scale_color_manual(values = c("#00BFC4", "#F8766D"), labels = c("Growth", "Decay")) +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"), labels = c("Growth", "Decay")) +
  theme_pubr(legend = c(0.1, 0.9), base_size = 13) ; FigX2

# Durga plot
d1 <- Durga::DurgaDiff(meandf, "meanforceperantN", "propmax")
p <- Durga::DurgaPlot(d1, contrasts = c("1 - 2", "2 - 1", "3 - 1", "4 - 1"), xlab = "", left.ylab = "Mean Force per ant (mN)", error.bars.type = "CI")

. Joining and leaving dynamics

Here we plot the joining and leaving dynamics for each replicate. Time is standardized so that each replicate starts at -1 and finishes at 1, where -1 to 0 indicates the growth phase, and 0 to 1 indicates the decay phase.

df1 <- df %>%
  group_by(replicate) %>%
  arrange(replicate, ptime) %>%
  mutate(
    changeNants = total.ants - lag(total.ants, 1) ## Find instances where number of ants changes
  )

df1$joinLeave <- ifelse(  ## Label each event as joining or leaving
  df1$changeNants == 1, "Join",
  ifelse(df1$changeNants == -1, "Leave", NA) 
)

df1$timeJoinLeave <- ifelse(
  !is.na(df1$joinLeave), df1$ptime, NA
)

JL_plot <- ggplot(df1, aes(ptime, forceN)) +
  geom_path() +
  geom_vline(data = df1[!is.na(df1$joinLeave),], aes(xintercept = timeJoinLeave, color = joinLeave), linetype = 2, na.rm = T) +
  facet_wrap(~replicate)

plotly::ggplotly(JL_plot) ## Interactive plot
050100050100050100-1.0-0.50.00.51.0050100-1.0-0.50.00.51.0-1.0-0.50.00.51.0-1.0-0.50.00.51.0
joinLeaveJoinLeaveptimeforceNCol210Col413_1Col413_2Col429Col517Col518Quel216Quel217Quel220_2Quel220_3Quel223Quel321_1Quel321_2Quel322Quel323

. Force change after joining or leaving event

Force change is evaluated over the 5 seconds following a joining or leaving event

timerange = 5 ## time of interest before/after each joining / leaving event 

df1 <- df1 %>%
  group_by(replicate, newarr) %>%
  mutate(
    duration = n()
  )

JL_events <- df1[!is.na(df1$joinLeave), ]

JL_forcechange = list()
for (row in 1:nrow(JL_events)) {

  event <- JL_events[row, ] 
  
  if (any(!is.na(df1[df1$replicate == event$replicate & 
                     between(df1$time, event$time - timerange, event$time - 1),]$joinLeave)) |
      any(!is.na(df1[df1$replicate == event$replicate & 
                     between(df1$time, event$time + 1, event$time + timerange),]$joinLeave)) |
      length(df1[df1$replicate == event$replicate & 
                     between(df1$time, event$time - timerange, event$time - 1),]$joinLeave) != timerange |
      length(df1[df1$replicate == event$replicate & 
                     between(df1$time, event$time + 1, event$time + timerange),]$joinLeave) != timerange) {
    next
    }
  
  force_pre <- mean(df1[df1$replicate == event$replicate & 
                     between(df1$time, event$time - timerange, event$time),]$forceN) 
  
  force_post <- (df1[df1$replicate == event$replicate & 
                     between(df1$time, event$time + 1, event$time + timerange),]$forceN)
  
  force_post <- force_post - force_pre[1]
  force_pre <- force_pre - force_pre[1]
  
  force <- c(force_pre, force_post)
  
  ants = event$total.ants
  ants = ifelse(between(ants, 1, 5), "1 - 5", ifelse(between(ants, 5, 10), "6 - 10", "10 +"))
  
  JL_forcechange[[length(JL_forcechange) + 1]] <- data.frame(
    replicate = rep(event$replicate[1], times = length(force_pre)),
    event = rep(event$joinLeave, times = length(force_pre)),
    ID = row,
    ants = ants,
    phase = rep(ifelse(event$ptime < 0, "Growth", "Decay"), times = length(force_pre)),
    force = force,
    time = seq(force) - 1
  )
}

JL_forcechange <- rbindlist(JL_forcechange)

JL_forcechange$phase <- factor(JL_forcechange$phase, levels = c("Growth", "Decay"))
JL_forcechange$ants <- factor(JL_forcechange$ants, levels = c("1 - 5", "6 - 10", "10 +"))

JL_forcechange_summ <- JL_forcechange %>%
  group_by(phase, event, time) %>%
  summarise(
    event = event[1],
    phase = phase[1],
    force_min = min(force),
    force_max = max(force),
    force_mean = mean(force),
    n = n(),
    SEM = sd(force)/sqrt(n()),
    sem_min = force_mean - SEM, 
    sem_max = force_mean + SEM,
    s = sd(force),
    margin = quantile(x = force, probs = 0.975)*sd(force)/sqrt(n()),
    ci_min = force_mean - margin, 
    ci_max = force_mean + margin
  )
## `summarise()` has grouped output by 'phase', 'event'. You can override using
## the `.groups` argument.
ForceTeam_JoinLeave <- ggplot(JL_forcechange_summ, aes(time, force_mean, fill = event)) + 
  geom_ribbon(aes(ymin = sem_min, ymax = sem_max), alpha = .4) + 
  geom_line(aes(color = event)) +
  geom_point(aes(color = event)) +
  scale_fill_viridis_d(option = "H", 
                       begin = 0.1, 
                       aesthetics = c('color', 'fill')) +
  labs(
    fill = "Event",
    color = "Event",
    x = "Time from event (s)",
    y = "Team force (mN)"
  ) +
  theme(
    legend.position = c(0.9, 0.85),
    legend.title = element_text(face = "plain")
  ) +
  facet_wrap(~phase); ForceTeam_JoinLeave
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

. Figure 2 - Main manuscript

## Note: Durga plots do not integrate well with ggplot and ggarrange. I save it as PNG and re-upload it to combine it with other plots.

# Create temp file
tmp_file <- tempfile(fileext = ".png")

# Open a PNG device and draw the Durga plot
png(filename = tmp_file, width = 1800, height = 1300, res = 300)
Durga::DurgaPlot(d1, contrasts = c("1 - 2" = "1 - 2", "2 - 1" = "2 - 1", "3 - 1" = "3 - 1", "4 - 1" = "4 - 1"), xlab = "", left.ylab = "Mean Force per ant (mN)", error.bars.type = "CI")
dev.off()
## png 
##   2
# Read the image and convert to grob
img <- readPNG(tmp_file)
p1 <- rasterGrob(img, interpolate = TRUE, 
                 width = unit(1, "npc"), height = unit(1.1, "npc"))

pplots <- list(FigX1, p1, FigX2, ForceTeam_JoinLeave)

# Combine plots
ggarrange(plotlist = pplots, nrow = 2, ncol = 2, labels = c("A", "C", "B", "D"))

. Bayesian model for force change following join/leave event

# Convert variables into factors
JL_forcechange$timeF <- as.factor(JL_forcechange$time)
JL_forcechange$ID <- as.factor(JL_forcechange$ID)
JL_forcechange$event <- as.factor(JL_forcechange$event)

JL_df <- JL_forcechange[JL_forcechange$time == 5, ]

# Data show extreme leptokurtosis and unimodal peak
ggplot(JL_df, aes(force)) + 
  geom_histogram() + 
  xlab("X")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# BRMS model
JL_df$force_z <- (JL_df$force - mean(JL_df$force))/sd(JL_df$force) # Centering and scaling

# Set robust prior
prior <- c(
  prior(student_t(3, 0, 10), class = "b"),  # fixed effects
  prior(gamma(2, 0.5), class = "nu")        # favors moderate tail thickness
)

# Model fit
model_brm <- brm(
  force_z ~ phase * event + (1 | replicate),
  data = JL_df,
  family = student(),
  prior = prior,
  iter = 4000, warmup = 1000, chains = 4, cores = 4,
  control = list(adapt_delta = 0.99)
)
## Compiling Stan program...
## Start sampling
# Model diagnostics
pp <- pp_check(model_brm, type = "dens_overlay")
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp + xlim(-20, 20)

pp_check(model_brm, type = "ecdf_overlay")         # Tails or peaks
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

pp_check(model_brm, type = "intervals")            # Predictive intervals per point
## Using all posterior draws for ppc type 'intervals' by default.

pp_check(model_brm, type = "stat", stat = "mad")   # Median absolute deviation
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

posterior_summary(model_brm, variable = "nu") # Estimated degrees of freedom
##    Estimate Est.Error     Q2.5    Q97.5
## nu 2.630123 0.7870871 1.531727 4.541379
mcmc_areas(as_draws_df(model_brm), pars = "nu") # Nu estimation looks robust for heavy tails

loo(model_brm)   # Pareto k diagnostics
## 
## Computed from 12000 by 184 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -216.6 15.3
## p_loo         9.2  0.5
## looic       433.2 30.7
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
### Results ###  
summary(model_brm)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: force_z ~ phase * event + (1 | replicate) 
##    Data: JL_df (Number of observations: 184) 
##   Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
##          total post-warmup draws = 12000
## 
## Multilevel Hyperparameters:
## ~replicate (Number of levels: 15) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.07      0.05     0.00     0.19 1.00     5667     5489
## 
## Regression Coefficients:
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                 0.55      0.12     0.32     0.79 1.00     6291
## phaseDecay               -0.51      0.17    -0.84    -0.19 1.00     7928
## eventLeave               -0.63      0.18    -1.00    -0.28 1.00     7849
## phaseDecay:eventLeave     0.27      0.22    -0.17     0.72 1.00     7355
##                       Tail_ESS
## Intercept                 8774
## phaseDecay                8359
## eventLeave                8705
## phaseDecay:eventLeave     7668
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.50      0.06     0.38     0.63 1.00     7984     7644
## nu        2.63      0.79     1.53     4.54 1.00     7619     7763
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## Hypothesis testing 

# Do joining events cause more force increase in the growth phase than in the decay phase?
hypothesis(model_brm, "phaseDecay < 0")
## Hypothesis Tests for class b:
##         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 (phaseDecay) < 0    -0.51      0.17    -0.79    -0.24    1089.91         1
##   Star
## 1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# Do leaving events cause more force decrease during decay phase than growth phase?
hypothesis(model_brm, "phaseDecay + phaseDecay:eventLeave < 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (phaseDecay+phase... < 0    -0.23      0.16     -0.5     0.02      14.13
##   Post.Prob Star
## 1      0.93     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# Does force increase for joining events during the decay phase?
hypothesis(model_brm, "Intercept + phaseDecay > 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Intercept+phaseD... > 0     0.04      0.12    -0.15     0.23        1.9
##   Post.Prob Star
## 1      0.65     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# Do leave events differ from 0 in the growth phase?
hypothesis(model_brm, "Intercept + eventLeave < 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Intercept+eventL... < 0    -0.08      0.15    -0.32     0.16       2.56
##   Post.Prob Star
## 1      0.72     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# Do leave events differ from 0 in the decay phase?
hypothesis(model_brm, "Intercept + phaseDecay + eventLeave + phaseDecay:eventLeave < 0")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Intercept+phaseD... < 0    -0.32      0.07    -0.44    -0.21        Inf
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

. Ant posture analysis

Visualize hindleg stretch based on ant position within the chain

allchainposdf$ant <- factor(allchainposdf$ant, levels = c('oneone','twoone','twotwo','threeone','threetwo','threethree'))

## Rename for ease of visualization
allchainposdf$nameplot <- ifelse(allchainposdf$ant == 'oneone', "Singletons", NA)
allchainposdf$nameplot <- ifelse(allchainposdf$ant == 'twoone', "2-ant-chain (front)", allchainposdf$nameplot)
allchainposdf$nameplot <- ifelse(allchainposdf$ant == 'twotwo', "2-ant-chain (back)", allchainposdf$nameplot)
allchainposdf$nameplot <- ifelse(allchainposdf$ant == 'threeone', "3-ant-chain (front)", allchainposdf$nameplot)
allchainposdf$nameplot <- ifelse(allchainposdf$ant == 'threetwo', "3-ant-chain (middle)", allchainposdf$nameplot)
allchainposdf$nameplot <- ifelse(allchainposdf$ant == 'threethree', "3-ant-chain (back)", allchainposdf$nameplot)

heat_df <- allchainposdf %>%
  group_by(chainsize, nameplot) %>%
  summarise(
    sd_stretch = sd(meanhindstretch, na.rm = TRUE), ## SD of hindleg stretch
    meanhindstretch = mean(meanhindstretch, na.rm = TRUE)) %>%  ## mean hindleg stretch
  ungroup()
## `summarise()` has grouped output by 'chainsize'. You can override using the
## `.groups` argument.
## Rearrange levels for plotting
heat_df$nameplot <- factor(heat_df$nameplot, levels = c(
  "Singletons",
  "2-ant-chain (front)", "2-ant-chain (back)",
  "3-ant-chain (front)", "3-ant-chain (middle)", "3-ant-chain (back)"
))

heat_df <- heat_df[order(heat_df$chainsize, heat_df$nameplot), ] 

# Add x/y positions to draw each ant in a spatial layout
heat_df$y <- as.numeric(as.factor(heat_df$chainsize))
heat_df$x <- ave(rep(1, nrow(heat_df)), heat_df$y, FUN = seq_along)

hindstretch_plot <- ggplot(heat_df, aes(x = x, y = -y)) +  # Flip y for top-down chain order
  geom_tile(aes(fill = meanhindstretch), width = 0.9, height = 0.9, alpha = 0.7) +
  geom_text(aes(label = paste0(round(meanhindstretch, 2), "\n± ", round(sd_stretch, 2))),
            size = 3.5, color = "black") +
  scale_fill_viridis(option = "B", name = "Extension (mm)", begin = 0.1) +
  scale_y_continuous(breaks = -unique(heat_df$y),
                     labels = paste0(unique(heat_df$chainsize), "-ant chain")) +
  theme(
    legend.position = c(0.7, 0.90),
    legend.title.position = "top",
    legend.direction = "horizontal",
    legend.key.width = unit(1.5, "cm"),      # widen color bar
    legend.key.height = unit(0.4, "cm"),     # adjust height
    legend.text = element_text(margin = margin(t = 5, b = 2)),  # space above/below labels
    legend.spacing.x = unit(0.5, "cm"),      # space between ticks
    legend.box.margin = margin(t = 5, r = 10, b = 5, l = 10),  # space around legend box
    plot.title = element_text(hjust = 0.5),
    axis.title.x = element_blank(),
    axis.title = element_blank(),
    text = element_text(size = 15)
  ) +
  scale_x_continuous(
    breaks = sort(unique(heat_df$x)),
    labels = c("Front", "Middle", "Rear")
  ) 

hindstretch_plot

. Figure 3 - Main manuscript

img <- readPNG("Images/ants.png")

g_img <- rasterGrob(img, interpolate = TRUE)

image_plot <- ggplot() + 
  annotation_custom(g_img, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf) +
  theme_void(); image_plot

ggarrange(image_plot, hindstretch_plot, ncol = 2, widths = c(0.8, 1), labels = c("A", "B"))

. Force of chains v. solo+chains

SoloVSoloChains_df <- meandf %>%
  group_by(replicate) %>%
  arrange(replicate, ptime) %>%
  mutate(
    phase = ifelse(ptime < 0, "Growth", "Decay")
  )

## Labelling
SoloVSoloChains_df$solos <- ifelse(  
  SoloVSoloChains_df$total.chains == 1, "Solos", 
  ifelse(SoloVSoloChains_df$chain1 > 0 & SoloVSoloChains_df$total.chains > SoloVSoloChains_df$chain1, "Solo + Chain", "Only chains"))

SoloVSoloChains_df$phase <- factor(SoloVSoloChains_df$phase, levels = c('Growth', 'Decay'))

SoloVSoloChains_df$solos <- factor(SoloVSoloChains_df$solos, levels = c("Solos", "Solo + Chain", "Only chains"))

SoloVSoloChains_df <- SoloVSoloChains_df[SoloVSoloChains_df$solos != "Solos",]

## Statistics

# Fit model
chainComb_mod <- glmmTMB(forceperantN ~ solos + total.ants + (1|replicate), SoloVSoloChains_df)

# Diagnostics
check_model(chainComb_mod) # Using performance
## `check_outliers()` does not yet support models of class `glmmTMB`.

res <- simulateResiduals(chainComb_mod) # Using DHARMa
testResiduals(res)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.049982, p-value = 0.3763
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.99848, p-value = 0.92
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 5, observations = 333, p-value = 0.1991
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.00489284 0.03469035
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                             0.01501502
# Results
summary(chainComb_mod)
##  Family: gaussian  ( identity )
## Formula:          forceperantN ~ solos + total.ants + (1 | replicate)
## Data: SoloVSoloChains_df
## 
##      AIC      BIC   logLik deviance df.resid 
##   1345.7   1364.8   -667.9   1335.7      328 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  replicate (Intercept) 1.523    1.234   
##  Residual              2.903    1.704   
## Number of obs: 333, groups:  replicate, 15
## 
## Dispersion estimate for gaussian family (sigma^2):  2.9 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       2.26683    0.42766   5.301 1.15e-07 ***
## solosOnly chains -0.28004    0.24550  -1.141    0.254    
## total.ants        0.22332    0.03176   7.031 2.05e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PLOTTING ## 
plot_model(chainComb_mod, "est") # Forest plot

SoloVSoloChains_plt <- ggplot(SoloVSoloChains_df[SoloVSoloChains_df$solos != "Solos",], aes(total.ants, forceperantN, color = solos)) + 
  geom_point(size = 3, alpha = 1) +
  facet_wrap(~phase) +
  labs(
    x = "Total number of ants",
    y = "Force per ant (mN)",
    color = "Arrangement"
  ) + 
  theme(
    text = element_text(size = 13),
    legend.position = "top"
  )

SoloVSoloChains_plt

. Number of ants over time, with Joining and Leaving events

## Use joining/leaving dataset
df2 <- df1 %>% 
  group_by(replicate) %>%
  arrange(replicate, ptime) 

df2$solos <- ifelse(  
  df2$total.ants == 1, "1", ifelse(between(df1$total.ants, 2, 5), "2 - 5", 
                                       ifelse(between(df2$total.ants, 6, 10), "6 - 10", 
                                              ifelse(between(df2$total.ants, 11, 13), "11 - 13", "13+"))))

## Order factor variable for plotting
df2$solos <- factor(df2$solos, levels = c("1", "2 - 5", "6 - 10", "11 - 13", "13+"))

# Preserve original names
unique_reps <- gtools::mixedsort(unique(df2$replicate)) 

df2 <- df2 %>%
  mutate(replicate = paste0("Replicate ", match(replicate, unique_reps)))

df2$replicate <- factor(df2$replicate, levels = unique(df2$replicate))

y_line <- min(df2$forceN, na.rm = TRUE) - 8
df_line <- data.frame(x = -1, xend = 1, y = y_line, yend = y_line)

# Plot
JL_time_plot <- ggplot() +
  geom_path(data = df2, aes(ptime, forceN, color = solos, group = replicate), linewidth = 0.8) +
  scale_color_brewer(palette = "Set1", name = "Nº of \nants") +
  
  new_scale_color() +
  
  geom_point(
    data = df2[!is.na(df2$joinLeave),],
    aes(x = ptime, y = -8, color = joinLeave),
    shape = "triangle",
    alpha = 0.7,
    size = 1.8
  ) +
  
  scale_color_manual(values = c("Join" = "black", "Leave" = "red"), name = "Event") +

  labs(x = "Normalised time", 
       y = "Team force (mN)", 
       color = "Event",
       shape = "Event") +
  
  ylim(-8, 135) +
  theme(text = element_text(size = 15)) +
  facet_wrap(~replicate)

JL_time_plot

. Force change between start and end of each arrangement

Density plot shows that, for each arrangement, force per ant increases during growth and decreases during decay. Only

df_forcechange_arrangement <- df %>%
  group_by(replicate, newarr) %>%
  mutate(
    duration = n(),
    timeFirst = mean(head(forceperantN), 3), ## Average force in the first 3 seconds of arrangement
    timeLast = mean(tail(forceperantN, 3)), ## Average force in the last 3 seconds of arrangement
    forcediff = timeLast - timeFirst,
    phase = ifelse(ptime < 0, "growth", "decay"),
    phase = factor(phase, levels = c("growth", "decay"))
  ) %>%
  slice_head()

force_diff_plot <- ggplot(df_forcechange_arrangement[df_forcechange_arrangement$duration > 10, ], aes(x = forcediff, color = phase)) +
  geom_density(linewidth = 1, key_glyph = draw_key_path) +
  labs(
    x = "Difference in force per ant between start and end of each arrangement (mN)",
    color = "Phase"
  ) +
  geom_vline(aes(xintercept = 0), linetype = 2, color = "black", linewidth = 0.5) +
  scale_color_manual(values = c("#00BFC4", "#F8766D"), labels = c("Growth", "Decay")) +
  theme_pubr(legend = c(0.1, 0.9),
             base_size = 13) +
  theme(strip.text = element_blank())

force_diff_plot

Thorax length v. force production

# Avoid name conflict
meandf <- meandf %>% select(-colony)

# Rename colonies
meandf$colony_short <- ifelse(startsWith(meandf$replicate, "Col"), substr(meandf$replicate, 1, 4), substr(meandf$replicate, 1, 5))

thoraxDF$colony_short <- str_remove(unique(thoraxDF$colony), "Thorax")

# Select relevant columns
thoraxDF <- thoraxDF[, c('colony_short', 'Length')]

thoraxDF <- thoraxDF[thoraxDF$colony_short %in% unique(meandf$colony_short), ]

# Thorax length is similar across colonies
ggplot(thoraxDF, aes(colony_short, Length, fill = colony_short)) +
  geom_boxplot()

# Merge datasets
thoraxDF <- thoraxDF %>%
  group_by(colony_short) %>%
  summarise(thorax_length = mean(Length))

meandf_wThorax <- merge(meandf, thoraxDF, by = "colony_short")

map = setNames(c("Colony 1", "Colony 2", "Colony 3", "Colony 4", "Colony 5"), c(unique(meandf_wThorax$colony_short)))

meandf_wThorax <- meandf_wThorax %>%
  mutate(colony_short = map[colony_short])


## PLOTTING ##
ggplot(meandf_wThorax, aes(thorax_length, meanforceperantN)) +
  geom_point(aes(color = colony_short)) +
  geom_smooth(method = "lm") +
  labs(
    color = "Colony", 
    y = "Mean force per ant (mN)",
    x = "Thorax length (cm)"
  )
## `geom_smooth()` using formula = 'y ~ x'

ggplot(meandf_wThorax, aes(thorax_length, chainsize)) +
  geom_point(aes(color = colony_short)) +
  geom_smooth(method = "lm") +
  labs(
    color = "Colony", 
    y = "Chain size (ants)",
    x = "Thorax length (cm)"
  )
## `geom_smooth()` using formula = 'y ~ x'

## STATISTICS ##
mod_thorax <- glmmTMB(meanforceperantN ~ thorax_length, meandf_wThorax)

# Diagnostics
res <- simulateResiduals(mod_thorax, n = 1000)
testResiduals(res)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.056451, p-value = 0.1964
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9988, p-value = 0.996
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 364, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.00000000 0.01008311
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0
# Results
summary(mod_thorax)
##  Family: gaussian  ( identity )
## Formula:          meanforceperantN ~ thorax_length
## Data: meandf_wThorax
## 
##      AIC      BIC   logLik deviance df.resid 
##   1593.6   1605.3   -793.8   1587.6      361 
## 
## 
## Dispersion estimate for gaussian family (sigma^2): 4.59 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)      2.753      7.928   0.347    0.728
## thorax_length    3.982     23.321   0.171    0.864